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Singular  perturbation  concepts  have  been  very  helpful  in  control 
and  systems  theory,  both  in  eliminating  difficulties  associated  with 
stiffness  and  in  reducing  dimensionality  (cf.  Kokotovic  et  a_l.  (1976) 
and  O'Malley  (1978)).  In  this  short  report,  we  wish  to  emphasize  the 
second  aspect,  especially  in  the  context  of  large-scale  problems.  Such 
an  approach  has  been  applied  to  power  system  modeling  (cf.  Kokotovic 
et  al.  (1979))  and  to  aircraft  maneuvering  and  structural  dynamics  (cf. 
Anderson  (1978)  and  Anderson  and  Hallauer  (1980)). 

Let  us  first  consider  initial  value  problems  on  a  bounded  interval, 
0  <  t  <  T,  for  linear  systems  of  the  form 
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x  =  Ax  +  Bu  ,  (1) 

where  A  and  B  are  constant  matrices  and  the  control  u(t)  is  considered, 
for  simplicity  in  the  present  discussion,  to  be  given. 

A  low-order  example  is  provided  by  a  model  of  the  longitudinal 
dynamics  of  an  F-8  aircraft  (cf.  Etkin  (1972)  and  Teneketzis  and 
Sandell  (1977)).  It  involves  four  states  x.  and  one  control  u,  the 
elevator  deflection.  Two  of  the  four  states,  velocity  variation  and 
flight  path  angle,  are  described  through  physical  intuition  as  being 
"primarily  (or  predominantly)  slow"  variables  compared  to  the  other 
"primarily  fast"  states,  angle  of  attack  and  pitch  rate.  It  is  a 
natural  temptation  in  such  a  situation  to  seek  to  approximate  the  solu¬ 
tion  away  from  t  =  0  by  neglecting  the  dynamics  of  the  primarily  fast 
T/ariables  and  to  integrate  only  the  resulting  initial  value  problem  for 
[the  slow  states.  (This  has  often  been  proposed  in  the  engineering 
literature,  cf.,  e.g.,  Harvey  and  Pope  (1976)).  The  trouble  with  this 
approach  is  that  the  physical  coordinates  are  not  actually  provided 
separately  as  slow  and  fast  modes,  in  contrast  to  the  situation  for  the 
traditional  singularly  perturbed  system  i  =  f(a,3,t),  e6  =  g(a,c,t) 
where  8  can  respond  on  a  much  faster  time-scale  than  a  as  £  *  0. 

The  analytical  solution  of  the  linear  system  (1)  is  commonly 
expressed  using  variation  of  parameters  and  the  matrix  exponential  so- 
lution  e  of  the  homogeneous  system.  Moler  and  van  Loan's  1978 
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article  "Nineteen  dubious  ways  to  compute  the  exponential  of  a  matrix" 
illustrates  the  difficulties  involved  in  computationally  exploiting 
specific  representations  for  this  fundamental  matrix.  The  difficulties 
are  compounded  when  the  dimensions  of  A  are  large.  We,  nonetheless, 

^  ifc 

recall  that  the  columns  of  e'  are  of  the  form  v. (t)e  where  v^  is  a 
vector  polynomial  in  t  in  the  span  of  the  eigenspace  corresponding  to 
the  eigenvalue  X^  of  A.  We  shall  seek  an  approximation  away  from  t  =  0 
which  avoids  doing  a  complete  eigenanalysis  of  A.  Our  approximation 
will  involve  neglecting  modes  corresponding  to  large  stable  eigenvalues 
of  A.  This  corresponds  to  the-  reduced  order  model  in  singular  pertur¬ 
bations  (cf.  O'Malley  (1974))  and  to  the  smooth  approximate  solution  of 
the  stiff  equations  literature  (cf.  Dahlquist  (1969)  and  Oden  (1971)). 

For  the  F-8  aircraft  model,  we  have  the  two  "slow"  (i.e.,  small  in 
magnitude)  eigenvalues  s.  ,  =  (-0.75  ±  i 7 . 6 )  *  10  2  compared  to  the  two 
"fast"  (i.e.,  large)  eigenvalues  f^  2  =  ~0*94  -  i3.0.  The  smallness  of 
the  parameter  u  =  [ s 2/ f 1 1  =  0*024  indicates  a  large  time-scale  separa¬ 
tion  between  slow  and  fast  modes,  and  the  smallness  of  e  =  -(1/Re  f^)/T 
=  1.06/T  for  a  long  enough  time  interval  0  <_  t  <_  T  provides  fast-mode 
stability  and  a  rapid  transition  to  pseudo-steady  state  behavior.  We 
shall  determine  a  second-order  nonstiff  model  for  the  solution  away 
from  t  =  0.  A  basic  question  is  what  initial  values  should  be  used  for 
the  second-order  problem.  (The  need  to  "project"  the  initial  vector 
occurs  in  many  related  situations  (cf.  O'Malley  and  Flaherty  (1980), 
Lentini  and  Keller  (1980),  and  de  Hoog  and  Weiss  (1979))).  Our 
approximation,  based  on  asymptotic  methods,  will  improve  as  the  para¬ 
meters  e  and  y  both  tend  toward  zero.  In  practice,  however,  there  is 
often  a  need  to  use  such  reduced-order  methods  even  when  these  para¬ 
meters  are  only  moderately  small.  In  contrast,  our  related  suggestion 
of  "equilibrating"  the  primarily  fast  states  will  produce  an  error 
which  persists  for  all  t  >  0  (cf.  Figure  1). 

In  the  applications  of  interest,  very  high  dimensional  systems  are 
common.  We  shall  seek  approximate  nonstiff  models  of  (differential) 
order  n^  <<  n,  the  dimension  of  x,  where  n^  is  the  number  of  "slow",  or 
small  in  magnitude,  eigenvalues  of  A  compared  to  its  other  (large) 
eigenvalues.  Systems,  for  which  the  eigenvalues  of  the  state  matrix  A 
can  be  so  divided,  will  be  called  two-time-scale .  Finer  subdivisions 
of  eigenvalues  into  slow,  fast,  and  very  fast  categories  could  be 
treated  by  repeating  our  procedure.  For  our  initial  value  problems, 
we'll  also  assume  fast-mode- stability,  i.e.  that  the  n£  =  n  -  n^  large 
eigenvalues  of  A  also  iave  large  negative  real  parts.  We  are  thereby 
eliminating  consideration  of  systems  with  large,  purely  imaginary 
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FIGURE  1.  F-8  aircraft  model:  Angle  of  attack 

(a  fast  variable)  vs.  time. 


eigenvalues.  Their  rapidly  oscillatory  solutions  require  techniques 
(cf.  Miranker  and  Wahba  (1976)  or  Petzold  (1978))  much  different  from 
the  boundary-layer  method  we  wish  to  develop. 

One  moderate-size  problem  we've  solved  is  a  recently  popular 
contest  problem  in  the  control  literature  (cf.  Sain  (1977)  and  De  Hoff 
and  Hall  (1979)).  It  is  a  sixteenth  order  model  for  a  turbofan  engine 
Indeed,  it  is  one  of  thirty  six  different  linearizations  used  to 
simulate  the  engine.  The  sixteen  eigenvalues  of  A  are  all  stable  and 
have  magnitudes  ranging  from  0.65  to  577.  Their  distribution  suggests 
using  either  a  third  or  a  fifth  order  reduced  model.  We  obtained  very 
good  approximations  in  both  cases  for  t  >  lOe,  with  e  being  determined 
in  terms  of  the  decay  rate  of  the  smallest  large  eigenvalue  (cf. 
O'Malley  and  Anderson  (1980)). 

To  proceed,  let  us  transform  the  constant  coefficient  problem 

x  =  Ax  +  Bu  (1 

into  a  new  problem 

y  =  Ay  +  Bu  ,  (2 


where  the  fast  and  slow  modes  of  the  homogeneous  problem  are  decoupled 


Specifically,  let 


y  -  Tx  , 


where 


a  =  Tat' 


B  =  7B 
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Here,  shall  have  the  small  eigenvalues  of  A  and  A^2  its  n£ 

(relatively)  large  eigenvalues  (with  large  negative  real  parts  as  well) 
We  note  that 

X  -K  ] 

1  |,  <5) 

!-L  I  +  LKj 
n2 


since,  e . g.  ,  T. 


f 1  0  1 

nl 

.  Partitioning  A  = 

A11  A12| 
A  A  1 

-L  I 

n2 

[  21  22j 

T  aT~^ 
lA'l 


A11  A12 

A11  ~  A12L 

A12 

.  °  A22. 

0 
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provided  the  n£  x  n^  matrix  L  satisfies  the  algebraic  Riccati  equation 


LA11  "  A22L  “  LA12L  +  A21  =  0  ■  (7) 

(The  connection  to  invariant  imbedding  is  signaled  by  the  entrance  of  a 
Riccati  equation.)  This  matrix  quadratic  equation  will  have  many 
solutions,  but  only  one  which  decouples  the  fast  modes  of  the  homo¬ 
geneous  problem,  as  in  (6)  .  Indeed,  if  the  spectrum  of  A^  coincides 

with  the  slow  eigenvalues  \  ,  i  =  l,...,ni,  of  A,  i.e.  if 
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it  is  easy  to  show  that 


where  the  n  *  matrix 


spans  the  n^  dimensional  slow  eigenspace 


of  A  (we  rearrange  the  rows  of  A,  if  necessary,  so  that  is 
invertible)  {cf.  Anderson  (1978),  Medanic  (1979),  and  van  Dooren  (1980) 
(which  develops  numerically  stable  computational  methods) ) .  Rather 
than  compute  this  eigenspace,  we  prefer  to  iterate  in  a  rearranged 
version  of  the  Riccati  equation  (7)  and  obtain  L  as  the  limit  of 
iterates 

Li+1  =  (A22  +  L1A12)-1<LiAu  *  A21)  .  (10) 

Not  unexpectedly,  this  procedure  relates  closely  to  Stewart  (1975)'s 
method  of  finding  the  dominant  eigenspace  of  A.  The  iteration  can  be 
shown  to  be  robust  with  respect  to  its  initialization,  and  its  rate  of 
convergence  is  proportional  to  the  time-scale  separation  parameter  y , 
i.e.  the  ratio  of  the  largest  slow  eigenvalue  in  magnitude  to  the 
smallest  fast  eigenvalue  in  magnitude. 

Using  the  L  so  obtained,  we  will  next  block  diagonalize  the  A  of 
(2)  provided  the  *  n ^  matrix  K  satisfies  the  linear  (Liapunov) 
equation 

KA22  "  AUK  +  A12  =  0  .  (11) 

Since  A^  and  A22  have  no  eigenvalues  in  common,  the  linear  system  has 
a  unique  solution  (cf.,  e.g..  Bellman  (1970))  which  can  be  found  as  the 
limit  of  the  iteration 

Vl  *  (XllKj  *  a12>*22  •  (12> 

Again,  the  procedure  has  a  rate  of  convergence  proportional  to  y. 

Having  now  determined  the  transforming  matrix  T  of  (4),  there 
remains  a  decoupled  system 

yx  =  A11y1  +  §xu  (13) 


y2  -  A22y2  +  B2u  , 


where  y  =  Tx  = 


for  y.  =  t.x  and  B.  =  r.B  for  ti  =  (I  +  KL  K)  and 

i  i  i  n  i 


and  t -  =  (LI  ).  We  note  that  this  decoupling  procedure  has  used  the 
n2 

the  two-time-scale  hypothesis,  but  not  the  hypothesis  requiring  fast¬ 
mode  stability. 


Since  we've  assumed  that  the  eigenvalues  of  A22  lie  far  into  the 
left  half  plane,  we  can  expect  y2  to  decay  rapidly  to  its  steady-state 


*2s  =  -X22S2U 


(15) 


(presuming  u  doesn't  vary  rapidly).  A  good  approximation  to  y  away 
from  t  *  0  will  then  follow  by  simply  integrating  the  nonstiff  system 
of  differential  equations  (3)  of  order  n^  with  initial  vector  y^CO)  =* 
t^x(O).  Returning  to  the  original  variables,  we  will  obtain  the  slow¬ 
mode  approximation 


(t)  =  r-i 


yl(t) 


2s 


(t) 


(16) 


It  is  generally  a  poor  approximation  near  t  =  0  because  the  "boundary 
layer  correction"  y2^  =  y2  ~  y2s  ne9lecte<l.  If  a  better  approxima¬ 
tion  is  needed  there,  one  can  integrate  the  linear  system 

y2f  =  ^22y2f  +  ^2s  with  =  T2X*°*  “  y2s^°^‘  Because  the  eigen¬ 

values  of  A22  are  large  and  (very)  stable,  we  will  need  to  use  a  short 
stepsize  only  on  a  short  initial  interval,  since  should  rapidly 

decay  to  zero.  This  would  fail  to  be  true  only  if  y2s  were  not  slowly- 
varying  (i.e.,  y2s  were  not  small)  for  t  >  0. 

A  more  substantial  problem  results  for  initial  value  problems  for 
the  time-varying  system 


x  =  A  ( t)  x  +  B  ( t)  u  ( t) 


(17) 


on  0  <  t  <  T.  Assuming  A  ( 0 )  is  two-time  scale  with  n-^  slow  modes,  we 
can  seek  a  transformation  T,  as  in  (4),  but  row  with  time-varying 
decoupling  matrices  L(t)  and  K(t)  determined  so  that  the  transformed 
problem  for 

y  =  f(t)x  (18) 


remains  two- time -scale  and  fast-mode  stable  with  the  same  fixed  number 
n^  of  (relatively)  small  eigenvalues  throughout  0  <_  t  <_  T. 

We  must  cautiously  note  that  eigenvalue  stability  does  not  imply 
stability  for  time-varying  systems  (cf . ,  e.g.,  Coppel  (1978)),  but  this 
is  nearly  true  for  certain  singularly  perturbed  systems  satisfying 
appropriate  stability  hypotheses  (cf.  the  statement  of  Tikhonov's 
theorem  in  Wasow  (1965)  or  Vasil 'eva  and  Butuzov  (1973)).  Kreiss 
(1978,  1979)  exhibits  counterexamples,  however,  for  systems  with  non¬ 
smooth  coefficients. 


Using  our  transformation  (4),  L(t)  will  block-triangularize  the 
system  matrix  A(t)  provided  it  now  satisfies  the  matrix  Riccati 
differential  equation 

L  *  -LAu(t)  +  A22(t)L  +  LA12(t)L  -  A2]_(t)  .  (19) 

If  L ( 0 )  =  0,  the  transformation  T^(0)  will  be  a  similarity  transforma¬ 
tion,  and  we  can  ask  (as  before)  that  the  eigenvalues  of  A^fO)  = 

Au(0)  -  A^2(0)L(0)  coincide  with  the  n^  slow  eigenvalues  of  A(0).  The 
initial  matrix  L(0)  can  then  be  obtained,  as  before,  by  iterating  in 
the  algebraic  Riccati  equation  L(0)  =  0.  We  will  then  obtain  L(t)  on 
0  <  t  <  T  by  integrating  the  resulting  initial  value  problem.  We  would 
stop  the  integration  and  abandon  our  procedure  if  we  encountered  a 
finite  escape  time,  or  if  the  transformed  problem  ceased  to  be  two- 
time-scale  and  fast-mode-stable  within  our  t  interval.  (Existence 
criteria  for  the  symmetric  matrix  Riccati  equation  commonly  encountered 
in  control  theory  are  known,  but  analogous  criteria  for  the  non-square 
problem  do  not  seem  to  be  available.)  The  possibility  for  intermittent 
reinitialization  should  be  investigated,  corresponding  to  the  reortho¬ 
normalizations  of  Scott  and  Watts  (1977). 

Knowing  an  appropriate  L(t),  the  matrix  K(t)  will  produce  a  block- 
diagonalization  of  the  state  matrix  provided  it  satisfies  the  linear 
differential  system 

K  =  Au(t)K  -  K  A  2  2  <  t )  -  A12(t)  .  (20) 

Since  fast-mode  stability  implies  that  the  eigenvalues  of  A22  = 

A22  +  *,al2  ^ave  ^ar<?e  negative  real  parts,  singular  perturbation 
concepts  (cf . ,  e.g.,  O'Malley  (.974))  suggest  that  the  solution  of  the 
initial  value  problem  for  K,  like  that  for  zz  =  z,  will  blow  up  for 
t  >  0.  Instead,  it  is  natural  to  integrate  a  terminal  value  problem 
for  K  and  expect  that  an  asymptotically  valid  approximation  will  result 
by  using  its  "pseudo-steady  state"  approximation,  i.e.  the  unique 
solution  of  the  algebraic  system 

K  ( t)  =  (Au(t)Kg(t)  -  A12(t)  )A22(t)  .  (21) 

The  need  for  nonuniform  convergence  at  terminal  time  will  be  eliminated 
by  picking  K(T)  =  0,  i.e.  K(T)  =  K  (T) . 

Since  the  homogeneous  variational  equation  for  L(t), 
i  *  -A^l  +  ;A22'  a^so  "singularly  perturbed",  but  opposite  in 
stability  to  the  system  for  K,  we  might  also  attempt  to  approximate 
L(t)  for  t  >  0  as  a  smooth  solution  L(t)  of  the  algebraic  Riccati 


equation  L(t)  =  0.  This  would  certainly  need  L(t)  to  be  slowly-varying 
with  the  eigenvalues  of  A22  +  LA^2  remaining  large,  and  far  into  the 
left  half-plane,  compared  to  those  of  A^  -  A^L.  Otherwise,  we  need 
to  actually  integrate  the  full  initial  value  problem  for  L(t). 

Now  note  that  the  fast-mode  stability  of  A^Ct)  implies  that  the 
slow-mode  or  pseudo-steady  state  approximation 

y2s(t)  =  “A^2 (t) B ( t) u (t)  (22) 

should  nicely  approximate  y2  for  t  >  0,  provided  y^s  remains  slowly- 
varying.  Thus,  we  are  finally  left  with  the  need  to  integrate  the  non¬ 
stiff  reduced-order  system  (13)  for  y^(t).  By  inverting  our  transfor¬ 
mation,  as  in  (16),  we  find  a  slow-mode  approximate  xs(t)  for  t  >  0. 

Note  that  all  initial  value  problems  for  x(t)  would  be  solved  in 
terms  of  an  n  x  n  fundamental  matrix  X(t)  for  the  homogeneous  system. 

In  the  two-time  scale  situation,  our  problem  splits  into  four  separate 
problems  for  the  n 2  x  n^  matrix  L(t),  for  the  n^  *  n2  matrix  K(t),  for 
the  n2  *  n^  fundamental  matrix  for  y 2 ( t ) ,  and  the  n^  *  n^  fundamental 
matrix  for  y^(t).  with  the  fast-mode  stability  assumption,  pseudo¬ 
steady  state  approximations  can  be  used  to  eliminate  the  differential 
systems  for  K,  Y2>  and  sometimes  L.  Thus,  a  substantial  order  reduction 
is  achieved  for  the  t  >  0  approximation. 

Substantially  more  detail  regarding  time-varying  problems  is 
contained  in  O'Malley  and  Anderson  (1980).  Various  related  problems 
might  be  treated  similarly.  If,  for  example,  a  two-time-scale  system 
produced  a  matrix  A22  with  only  moderate-sized  eigenvalues,  one  might 
consider  y^  to  be  approximately  constant  on  finite  T  intervals,  so  that 
only  a  differential  system  for  y2  need  be  integrated.  Likewise,  if  the 
eigenvalues  of  A22  are  large  with  both  large  negative  and  large  posi¬ 
tive  real  parts,  one  might  seek  reduced-order  approximations  for 
appropriate  two-point  problems.  Finally,  extensions  of  these  ideas  to 
nonlinear  problems  must  be  sought. 
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